(* Mathematica Package *)
(*Clear["TransferMatrixFormalism`*"];*)
BeginPackage["TransferMatrixFormalism`"];

Unprotect["`*"]; 
ClearAll["`*"];

RefractionAngle::usage="will be imported";
FresnelA::usage="will be imported";
RefractionMv::usage="will be imported";
PhaseMv::usage="will be imported";
AbelesS::usage="will be imported";
AbelesSIv::usage="will be imported";
AbelesSIIv::usage="will be imported";
TransmitTopS::usage="will be imported";
ReflectTopS::usage="will be imported";
TransmitBottomS::usage="will be imported";
ReflectBottomS::usage="will be imported";
InternalTransferC::usage="will be imported";
ExternalTransferC::usage="will be imported";
(*The "will be imported" method allows Wolfram Workbench to recognize that
these symbols are in the main context (not `Private`) while keeping
symbol::usage mess in a separate file.  Then, symbol documentation can be
autogenerated using Usages defined there:*)

Get[FileNameJoin[{"TransferMatrixFormalism","Usages", "TransferMatrixFormalism.m"}]]


Begin["`Private`"]

Unprotect["`*"]; 
ClearAll["`*"];


(*Snell's Law*)
RefractionAngle = ArcSin[ #2 / #3 Sin[#1]] &;(*[\[Theta]1_,n1_,n2_]*)


(*Optical Admittance*)
ys = -#2 Cos[#1] &;

yp = #2/Cos[#1] &;

cp = Cos[#1]/Cos[#2] &;


(*Internally Implemented Fresnel amplitude coefficients for reflection and transmission*)
rs = Evaluate[(ys[#1, #3] - ys[#2, #4]) / (ys[#1, #3] + ys[#2, #4])]&;

rp = Evaluate[(yp[#1, #3] - yp[#2, #4]) / (yp[#1, #3] + yp[#2, #4])]&;
	
ts = Evaluate[(2 * ys[#1, #3]) / (ys[#1, #3] + ys[#2, #4])]&;

tp = Evaluate[cp[#1, #2] * (2 * yp[#1, #3]) / (yp[#1, #3] + yp[#2, #4])]&;

(*End user Fresnel Amplitude Coefficients*)
headQ = MemberQ[{"rs","rp","ts","tp"},#]&

FresnelA[head_?headQ, thetai_, ni_, nj_] := 
	ToExpression["TransferMatrixFormalism`Private`"<>head][thetai, RefractionAngle[thetai,ni, nj], ni, nj];
	
FresnelA[head_?headQ, thetai_, thetaj_, ni_, nj_] := 
	ToExpression["TransferMatrixFormalism`Private`"<>head][thetai, thetaj, ni, nj];


(*End-user functions for calculating transmission and reflection through S, SI, or SIIv*)
SDimQ = (Dimensions[#] == {2,2})&;
TransmitTopS[s_?SDimQ] := 1/s[[1, 1]];
ReflectTopS[s_?SDimQ] := s[[2, 1]]/s[[1, 1]];
TransmitBottomS[s_?SDimQ] := Det[s]/s[[1, 1]];
ReflectBottomS[s_?SDimQ] := -(s[[1, 2]]/s[[1, 1]]);

(*Internally implemented functions for calculating transmission and reflection through S, SI, or SIIv*)
tTOP[s_] := 1/s[[1, 1]];
rTOP[s_] := s[[2, 1]]/s[[1, 1]];
tBOT[s_] := Det[s]/s[[1, 1]];
rBOT[s_] := -(s[[1, 2]]/s[[1, 1]]);




(*Matrix of Refraction*)
Get[FileNameJoin[{"TransferMatrixFormalism","DownValues", "RefractionMv.m"}]]

(*Phase Matrix*)
Get[FileNameJoin[{"TransferMatrixFormalism","DownValues", "PhaseMv.m"}]];


PolQ=MemberQ[{"s","p"},#]&;
(*Total System Transfer Matrix*)
Get[FileNameJoin[{"TransferMatrixFormalism","DownValues", "AbelesS.m"}]]

(*Partial System Transfer Matrix: Subsystem I*)
Get[FileNameJoin[{"TransferMatrixFormalism","DownValues", "AbelesSIv.m"}]]

(*Partial System Transfer Matrix: Subsystem II*)
Get[FileNameJoin[{"TransferMatrixFormalism","DownValues", "AbelesSIIv.m"}]]


denomv[AbelesSIv_, PhaseMv_, AbelesSIIv_] := (1 - rBOT[AbelesSIv] rTOP[AbelesSIIv] (PhaseMv[[2, 2]])^2);
(*Internal Transfer Coefficients*)
Get[FileNameJoin[{"TransferMatrixFormalism","DownValues", "InternalTransferC.m"}]]

(*External Transfer Coefficient*)
Get[FileNameJoin[{"TransferMatrixFormalism","DownValues", "ExternalTransferC.m"}]]


Protect["`*"]; 
End[]

Protect["`*"]; 
EndPackage[]